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Materials characterized by spatially homogeneous elastic moduli undergo affine distortions when 
subjected to external stress at their boundaries, i.e., their displacements u(x) from a uniform refer- 
ence state grow linearly with position x, and their strains are spatially constant. Many materials, 
including all macroscopically isotropic amorphous ones, have elastic moduli that vary randomly 
with position, and they necessarily undergo nonaffine distortions in response to external stress. We 
study general aspects of nonaffine response and correlation using analytic calculations and numeri- 
cal simulations. We define nonaffine displacements u'(x) as the difference between u(x) and affine 
displacements, and we investigate the nonaffinity correlation function Q(x.) = {[u'(x) — u'(0)] 2 } and 
related functions. We introduce four model random systems with random elastic moduli induced 
by locally random spring constants (none of which are infinite), by random coordination number, 
by random stress, or by any combination of these. We show analytically and numerically that 5(x) 
scales as Alxl - ^ -2 ' where the amplitude A is proportional to the variance of local elastic moduli 
regardless of the origin of their randomness. We show that the driving force for nonaffine displace- 
ments is a spatial derivative of the random elastic constant tensor times the constant affine strain. 
Random stress by itself does not drive nonaffine response, though the randomness in elastic moduli 
it may generate does. We study models with both short and long-range correlations in random 
elastic moduli. 

PACS numbers: 62.20.Dc, 83.50.-v, 83.80.Fg 



I. INTRODUCTION 

In the classical theory of elasticity Q, 0, 0, , an elas- 
tic material is viewed as a spatially homogeneous medium 
characterized by a spatially constant elastic-modulus ten- 
sor Kijki- When such a medium is subjected to uniform 
stresses at its boundaries, it will undergo a homogeneous 
deformation with a constant strain. Such homogeneous 
deformations are called affine. This picture of affine 
strain is generally valid at length scales large compared 
to any characteristic inhomogeneities: displacements av- 
eraged over a sufficiently large volume are affine (at least 
in dimensions greater than two). It applies not only to 
regular periodic crystals, but also to polycrystalline ma- 
terials like a typical bar of steel. At more microscopic 
scales, however, individual particles in an elastic medium 
do not necessarily follow trajectories defined by uniform 
strain in response to external stress: they undergo non- 
affine rather than affine displacements. The only systems 
that are guaranteed to exhibit affine distortions at the 
microscopic scale are periodic solids with a single atom 
per unit cell. Atoms within a multi-atom unit cell of 
a periodic solid will in general undergo nonaffine distor- 
tions p| , and atoms in random and amorphous solids will 
certainly undergo nonaffine distortions. Such distortions 
can lead to substantial corrections to the Born-Huang Q 
expression for macroscopic elastic moduli. 

Research on fragile ® R, 11 | oL Eo |. g ranular [jTlH^. 
crosslinked polymeric IT3i. fl5l |l6L Tl7L N~8| . and biolog- 
ical materials fl9l I20L I2lll22l |23|. particularly in small 
samples, has sparked a renewed interest in the nature of 



nonaffine response and its ramifications. Liu and Langer 
Q introduced various measures of nonaffinity, in partic- 
ular the mean-square deviation from affinity of individ- 
ual particles in model foams subjected to shear. Tan- 
guy et al. in their simulation of amorphous sys- 
tems of Lennard- Jones beads found substantial nonaffine 
response and a resultant size-dependence to the macro- 
scopic elastic moduli. Lemaitre and Maloney [25( relate 
nonaffinity to a random force field induced by an initial 
affine response. Head et al. [20l l2ll studied models 
of crosslinked semi-flexible rods in two-dimensions and 
found two types of behavior depending on the density of 
rods. In dense systems, the response is close to affine 
and is dominated by rod compression, whereas in more 
dilute systems, the response is strongly nonaffine and 
dominated by rod bending. 

The recent work discussed above provides valuable in- 
sight into the nature of nonaffine response. It does not, 
however, provide a general framework in which to de- 
scribe it. In this paper, we provide a such a frame- 
work for describing the long-wavelength properties of 
nonaffinity, and we verify its validity with numerical 
calculations of these properties on a number of zero- 
temperature central-force lattice models specifically de- 
signed to demonstrate our ideas. Our hope is that this 
framework will prove a useful tool for studying more re- 
alistic models of amorphous glasses, granular material, 
and jammed systems, particularly at zero temperature 
just above the jamming transition [26L 1211 128|. We are 
currently applying them to jamme d sy stems |29j | and to 
networks of semi-flexible polymers [30|. 
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Though nonaffinity concerns the displacement of indi- 
vidual particles at the microscopic scale, we show that 
general aspects of nonaffine response in random and 
amorphous systems can be described in terms of a con- 
tinuum elastic model characterized by a local elastic- 
modulus tensor Kijki(x) at point x, consisting of a spa- 
tially uniform average part Kijki and a locally fluctu- 
ating part SKijki(x), and possibly a local stress tensor 
&ij (x) with vanishing mean. We show that under stress 
leading to a macroscopic strain 7^, the random part of 
the elastic-modulus tensor, in conjunction with the strain 
7y, acts as a source of nonaffine displacement u-(x) pro- 
portional to djSKijkiJkl- F° r small SKijki and 7^, the 
Fourier transform of the correlation function dj (x, 0) of 
the displacement u'(x) can be expressed schematically 
as 7 2 A K (q) / (q 2 K 2 ) where A K (q) represents the Fourier 
transform of relevant components of the correlation func- 
tion of the random part of the elastic-modulus tensor 
and K represents the average elastic-modulus tensor. At 
length scales large compared to the correlation length 
£ of the random elastic modulus, A^(q) is a constant 
A^, and the nonaffinity correlation function in d dimen- 
sions scales as (A K / K 2 )~/ 2 \x\~( d ~ 2 \ which exhibits, in 
particular, a logarithmic divergence in two dimensions; 
at length scales smaller than £, A K (q) ~ q~^, where 
<fi can be viewed as a critical exponent, and the non- 
affinity correlation function scales as |x|^+ 2 ~ d for tj> < d. 
For simplicity, we focus on zero-temperature systems. 
Our analytic approach is, however, easily generalized to 
nonzero temperature in systems with unbreakable bonds. 
At nonzero temperature, the dominant, long-distance be- 
havior of nonaffinity correlation functions is the same as 
at zero temperature. 

Our numerical studies were carried out on systems 
composed of sites either on regular periodic lattices or 
on random lattices constructed by sampling a Lennard- 
Jones liquid and connecting nearest-neighbor sites with 
unbreakable central force springs. We allowed the spring 
constants of the springs, their preferred lengths, or both 
to vary randomly. The local elastic modulus at a par- 
ticular site in these models depends on the strength and 
length of springs connected to that site as well as on the 
number of springs connected to it. Thus, a periodic lat- 
tice with random spring constants and an amorphous lat- 
tice with random site coordination numbers both have a 
random local elastic constant. Their nonaffinity correla- 
tion function should, therefore, exhibit similar behavior, 
as our calculations and simulations verify. It is important 
to note that macroscopically isotropic systems are always 
amorphous and, therefore, always have a random elastic- 
modulus tensor and exhibit nonaffine response. For sim- 
plicity, we do not consider systems in which any spring 
is infinitely rigid (i.e., has an infinite spring constant). 
With appropriate coarse graining of 6Kijki{x), however, 
our primary analytical results are expected to apply to 
this more general case. 

The outline of this paper is as follows. In Sec. [n] we 
derive familiar formulae for the elastic energy of central- 



force lattices and introduce our continuum model, giving 
special attention to the nature of random stresses. In 
Sec. IIIII we use the continuum model to calculate non- 
affine response functions in different dimensions for sys- 
tems with random elastic moduli with both short- and 
long-range correlations and with random stress tensors 
relative to a uniform state, and we calculate the corre- 
lation function of local rotations induced by nonaffine 
distortions. In Sec. IIVI we present numerical results for 
the four model systems we consider: periodic lattices 
with random elastic constants without (Model A) and 
with (Model B) random stress, and amorphous lattices 
with random elastic constants without (Model C) and 
with (Model D) random stress. Four appendices present 
calculational details: Appendix ^ derives the indepen- 
dent components of the 8th rank modulus correlator in 
an isotropic medium, Add. IbI calculates the general form 
of the nonaffinity correlation function as a function of 
wavevector, App. [21 calculates the asymptotic forms as 
a function of separation x of the nonaffinity correlation 
function, and App. IU1 calculates the correlation function 
of local vorticity. 

II. MODELS AND DEFINITIONS 

A. Notation and Model Energy 

We consider model elastic networks in which particles 
occupy sites on periodic or random lattices in their force- 
free equilibrium state. Thus, particle £ is at lattice posi- 
tion R^o in equilibrium. When the lattices are distorted, 
particle £ undergoes a displacement to a new position 

R^ = R ffl + Uf. (2.1) 

We will refer to the equilibrium lattice, with lattice posi- 
tions R^o, as the reference lattice or reference space, 
and the space into which the lattice is distorted via the 
displacements as the target space. Pairs of parti- 
cles £ and £' are connected by unbreakable central-force 
springs on the bond b =< £',£ >. The coordination 
number of each particle (or site) is equal to the number 
of particles (or sites) to which it is connected by bonds. 
The potential energy, Vb(Rb), of the spring on bond b 
depends only on the magnitude, 

Rb = |Rf — Rf|, (2-2) 

of the vector connecting particles £ and £'. The total 
potential energy is thus 

U T = ^2v b (R b ) = i^I/<^>(|R^ -R,|). (2.3) 

b 1,1' 

We will consider anharmonic potentials 

V b = ^h{R b - R bR ) 2 + -g b (R b - R bR )\ (2.4) 
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with both harmonic and quartic components, where RbR 
is the rest length of bond b. We assume that both kb and 
gb are finite. The harmonic limit is obtained when the 
quartic coefficient gb vanishes, in which case, kb is the 
harmonic spring constant. 

We will only study systems in which there is an equi- 
librium reference state with particle positions {R^o} in 
which the force on each site is zero. The length Rbo = 
\RfO — Rio I °f each bond b in this configuration does not 
have to coincide with its rest length RbR- As we shall 
see in more detail shortly, it is possible to have the total 
force on every site be zero but still have nonzero forces 
on each bond. 

The potential energy of the lattice can be expanded in 
terms of the discrete lattice nonlinear strain [3, 

«6 = ~ (Rl - Rbo) = R w • Au & + \ ( Au 6 ' Au b) ( 2 -5) 

relative to the reference state, where Au b = — Up. 
The discrete strain variable, v b , is by construction in- 
variant with respect to rigid rotations of the sample, i.e., 



it is invariant under Rp 



UijRej, where Uij is any 



independent rotation matrix. To second order in Vb in 
an expansion about a reference lattice with lattice sites 
R^o, the potential energy is 



b 



(2.6) 



where F(b) = |F(6)| is the magnitude of the force, 

F(6) = -V b \R b o)Rb /R b o, (2.7) 
acting on bond b and 



k(b) = V b ^(R bo )-R b - 1 V b \ ) (R b o) 



(2- 



is the effective spring constant of bond 6, which reduces 
to kb when Rbo — RbR- k(b) is never infinite because we 
we assume kb and gb are finite. The equilibrium bond- 
length Rbo for each bond is determined by the condition 
that the total force at each site i vanish at ui = 0: 



dAUj 



due 



u f =0 



= e',e>) 



(2.9) 



where e b Qi = Rboi/Rbo is the unit vector directed along 
bond b. Thus the harmonic potential on each bond de- 
composes into a parallel part, proportional to V b , di- 
rected along the bond and a transverse part, propor- 
tional to R b Q V b , directed perpendicular to the bond. The 
transverse part vanishes when the force on the bond van- 
ishes. 

The harmonic energy A/7^ ar does not preserve the in- 
variance with respect to arbitrary rotations of the full 
nonlinear strain energy AUt of Eq. I|2.6|) . under which 



Auu — > Ai 



{U lj - 5 i3 )R b0j + UijAubj, (2.11) 



where Uu is a rotation matrix. It does, however preserve 
this invariance up to order 9 2 but not order 9 2 Au b and 
9(Au b ) 2 , where 9 is a rotation angle. For small 0, 



Au' b = Au b + 0x R b0 + O{0 2 ,9Au b ), 



(2.12) 



and e h0 • Au' b = e b0 ■ Au b + 0(9 2 ,9Au b ). Thus, the 
part of the harmonic energy arising from the k(b) term 
in Eq. I|2.6|l is invariant to the order stated above. The 
invariance of the force term of Eq. I|2.6|l is more subtle. 
Under the above transformation of Eq. (|2.12|l . (Au' b ) 2 = 
(Au b ) 2 + 26 x R b ■ Au b + (0 x R b ) 2 + 0{6 2 Au b , e{Au h ) 2 ), 
and it would seem that there are terms of order 9, and 9 2 
in A£/^. al . These terms vanish, however, upon summation 
over £ and £' because of the equilibrium force condition 
of Eq. 12.9[> . Thus, the full AU^ is invariant under 
rotations up to order 9 2 . 



B. Definition of Models 

We will consider the following simple models of 
random lattices. 

Model A: Random, zero-force bonds on a 
periodic lattice. In this model, all sites lie on a 
periodic Bravais lattice with all bond lengths constant 
and equal to Rbo, and the rest length RbR of each bond 
is equal to Rbo- The force F(b) on each bond is zero, 
but the spring constant fc b and other properties of the 
potential V b can vary from site to site. Each lattice site 
has the same coordination number. 



This equilibrium condition only requires that the total 
force on each site, arising from all of the springs attached 
to it, be zero [3l|]. It does not require that the force F(b) 
be equal to zero on every bond b. 

In equilibrium, when Eq. I|2.9|l is satisfied, the part of 
v b linear in Au b disappears from AUt- In this case, it is 
customary to express AUt to harmonic order in Au b : 
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Model B: Random, finite-force bonds on an 
originally periodic lattice. In this model, sites are 
originally on a regular periodic lattice, but rest bond 
lengths RbR are not equal to the initial constant bond 
length on the lattice. Sites in this model will relax to 
positions Rio with bond lengths Rbo = |Ri'o — Rfol such 
that the force F(^) at each site I is zero but the force 
F(6) exerted by each bond b is in general not. This 
model has random stresses and, as we shall see, random 
elastic moduli as well. The bond vectors R b o and spring 
eboiebOj)]AubiAubj, constant fc b are random variables, but the coordination 
number of each site is not. Random stresses in an 
(2.10) originally periodic lattice necessarily induce randomness 
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in the elastic moduli relative to the relaxed lattices with 
zero force at each site. 

Model C: Random, zero-force bonds on a 
random lattice. In this model, lattice sites are 
at random positions and have random coordination 
numbers. The equilibrium length Rbo varies from bond 
to bond. The rest length RbR of each bond is equal 
to its equilibrium length so that the force F(b) of each 
bond is zero. This model, which is meant to describe 
an amorphous material, is macroscopically but not 
microscopically homogeneous and isotropic. 

Model D: Random finite-force bonds on a 
random lattice. This is the most general model, and it 
is the one that provides the best description of glassy and 
random granular materials. In it, the rest lengths RbR, 
the spring constants kb, and the coordination number 
are all random variables. Like Model C, this model 
describes macroscopically isotropic and homogeneous 
amorphous material. 



Though Models A, B, and C can be viewed as sub- 
sets of the most general model D, we find it useful to 
treat them as distinct models because they each isolate 
separate causes of randomness in the local elastic mod- 
ulus or stress. One of our goals, for example, is to show 
analytically and numerically that the non-affinity corre- 
lations arising from structural randomness in models C 
and D have exactly the same form as those arising from 
the more controlled periodic models A and B. Another 
is to study the different effects of random elastic moduli 
and random stress. 

In all of these models the random elastic-modulus ten- 
sor can in principle exhibit either short- or long-range 
correlations in space. To investigate the effects of such 
long-range correlations, we explicitly construct spring 
constant distributions with long-range correlations in 
model A. We will also find evidence of long-range correla- 
tions in model C when the reference lattice has correlated 
crystalline domains. 



C. Continuum Models 

In the continuum limit, when spatial variations are 
slow on a scale set by the lattice spacing, the equilib- 
rium lattice positions become continuous positions x in 
the reference space: R^o — * x; and the target-space po- 
sition and displacement vectors become functions of x: 
R^ — » R(x) and — > u(x). In this limit, the lattice 
strain Vb becomes 



where 



u ij( x ) = n(diUj + djU t + diU ■ dju) 



(2.13) 



(2.14) 



is the full Green-Saint Venant Lagrangian nonlinear 
strain Q, 0, 0] , which is invariant with respect to rigid 
rotations in the target space [i.e., with respect to rigid 
rotations of R(x)]. Sums over lattice sites of the form 
^2lS(£), for any function S(£), can be replaced by in- 
tegrals J d d xS(x)/v(x) where u(x) is the volume of the 
Voronoi cell centered at position x = R^o- The contin- 
uum energy is then 



H= d d x 



where 



-K ijM (x)My (x)u kt (x) + aij (x)?iy (x) 



M*) = -2^x)E^ 



(2.15) 

boj\b=<e>,i> (2-16) 



is a local symmetric stress tensor at x where the sum 
over t 1 is over all bonds with one end at I and 



Kijki (x) 



1 



2«(x) 



^ k(b)R b Q Rb0iRb0jRb0kRb0l\b=<£' ,£> 



(2-17) 

is the local elastic-modulus tensor [33. Because it de- 
pends only on the full nonlinear strain Ujj(x), the con- 
tinuum energy TL of Eq. (|2.15|l is invariant with respect to 
rigid rotations in the target space. This is a direct result 
of the fact that we consider only internal forces between 
particles. The stress tensor cr^ (x) is generated by these 
internal forces, and as a result, it multiplies «y in TL. It is 
necessarily symmetric, and it transforms like a tensor in 
the reference space. (It is not, however, the second Piola- 
Kirchoff tensor [3j|, a{J — 5TL/5uij(x) — KijkiUki +5^, 
which also transforms in this way.) External stresses, on 
the other hand, specify a force direction in the target 
space and couple to the linear part of the strain. 

Since Kij k t(x) in Eq. (|2.17(l arises from central forces 
on bonds, it and its average over randomness obey the 
Cauchy relations 0,0, Ki jk i(x) = K ik ji(x) = K ajk {x), 
in addition to the more general symmetry relations, 
2£yfcj(x) = K jik i(x) = Kijiu(x) = K klij (x). The Cauchy 
relations reduce the number of independent elastic mod- 
uli in the average modulus Kij k i = (l£yjy(x)) below 
the maximum number permitted for a given point-group 
symmetry (for the lowest symmetry, from 21 to 15). In 
particular, they reduce the number of independent mod- 
uli in isotropic and hexagonal systems from two to one, 
setting the Lame coefficients A and /x equal to each other. 
In our analytical calculations, we will, however, treat A 
and fj, as independent. The Cauchy limit is easily ob- 
tained by setting A = fx. 

The stress tensor (x) is generated entirely by inter- 
nal forces on bonds. The elastic-modulus tensor Kij k i(yi) 
depends on the local effective spring constant k(b), the 
length and direction of the bond vectors R&o, and the site 
coordination number; and it will be a random function of 
position if any of these variables are random functions of 
position. Thus Kij k i(x) is a random function of position 
in Models A to D. The stress tensor &ij(x) is nonzero 
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only if the bond forces are nonzero. It is thus a random 
function of position only in Models B and D. 

We require that the continuum limit of our lattice mod- 
els be in mechanical equilibrium when u(x) = 0. This 
means that the linear variation of TL with respect to u(x) 
must be zero, i.e., that 

SH = J d d xa ij (x)d j Su i (x) = (2.18) 

for any 5ui(x). 5ui(x) can be decomposed into a con- 
stant strain part and a part whose average strain van- 
ishes: 6ui(x) = S-fijXj + Ju'(x) where f d d djSu' i (x) = 0. 
Equilibrium with respect to variations in 7.^ implies that 
the spatial average of <jy is zero. Equilibrium with re- 
spect to <5u'(x) implies that when x is in the interior of 
the sample, 

/<(x) = djaij (x) = djOji = 0, (2.19) 

where f is the force density that is a vector in the target 
space. In addition, J dSj&ij (x)du^ (x) = for any Suj(x), 
where the integral is over the surface of the sample, im- 
plying that ^ (xb) = for points on the surface. 

Thus, we see that equilibrium conditions in the refer- 
ence space impose stringent constraints on the random 
stress tensor 5y(x): its spatial average must be zero, its 
values on sample surfaces must be zero, and it must be 
purely transverse, i.e., it must have no longitudinal com- 
ponents parallel to the gradient operator. Though the 
linear part of u y does not contribute to the stress term 
in 7i, the nonlinear part still does, and 7i can be written 
as 

H = i J d d x[K ijk i(x)u ij (x)u k i(x) 

+aij(x)d i Uk(x)djUk( f x.)]. (2.20) 

Because of the constraints on cr^ , this free energy is iden- 
tical to that of Eq. (|2.15() . It is invariant with respect to 
rotations in the target space even though it is written so 
that the explicit dependence on the rotationally invariant 
strain is not so evident |34|. 

As we have seen, the spatial average of &ij (x) is zero; it 
only has a random fluctuating part in models we consider. 
The elastic- modulus tensor Kij k i(x), on the other hand, 
has an average part and a random part with zero mean: 

K ijkl (x) = K ijkl + 5K ijk i(x). (2.21) 

We will view both a%j (x) and SKijM (x) as quenched ran- 
dom variables with zero mean. 



III. STRAINS AND NONAFFINITY 

Consider a reference elastic body in the shape of a 
regular parallelepiped. When such a body is subjected 
to stresses that are uniform across each of its faces, it 




(a) (b) 

FIG. 1: Sheared elastic medium with nonamne displacements, 
(a) Unsheared reference state, (b) Sheared state with non- 
amne displacements. Under affine distortions, points on on 
the vertical dotted lines in (a) would map to points on the 
slanted dotted lines parallel to the left and right of bound- 
aries of the sheared sample in (b); under nonafnne distortion, 
they do not. [Color online] 

will undergo a strain deformation in which its boundary 
sites at positions xg distort to new positions 

Ri{x B ) = AijXBj, (3.1) 

where A y - is the deformation gradient tensor [33| . If the 
medium is spatially homogeneous, then A y - = Sij + 7^ 
determines the displacements of all points in the medium: 
Ri(x) = AijXj or itj(x) = jijXj. Such a distortion is 
called affine. In inhomogeneous elastic media, there will 
be local deviations from affinity [Fig. ^ described by a 
displacement variable u'(x) defined via 

Ri(x) = AijXj +u' i (x) (3.2) 

or, equivalently, 

iij(x) = jijXj + u-(x) (3.3) 
Uijix) w 7?.£j + (diu'j + dju't 

+li P dju' p + i jp diu' p )/2, (3.4) 

where the final equation contains only terms up to linear 
order in u' and where 7^ = (7^ + "fji + "fik'Jjk)/'^- Since 
distortions at the boundary are constrained to satisfy Eq. 
(|3.1() . !ij(xfl) is zero for all points xgon the boundary. It 
is often useful to consider periodic boundary conditions 
in which u'(x) has the same value (possibly not zero) 
on opposite sides of the parallelepiped. This condition 
implies 

J d d xd ] u' i {x) = jdSju'i = 0. (3.5) 



A. Nonaffinity in Id 

To develop quantitative measures of nonaffinity, it is 
useful to consider a simple one-dimensional model, which 
can be solved exactly. We study a one-dimensional pe- 
riodic lattice, depicted in Fig. |3 with sites labelled by 
I = 0, ■■■,N, whose equilibrium positions are Rio = ai, 
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FIG. 2: Schematic diagram of nonaffine distortion in a one- 
dimensional lattice with random spring constants. The top 
figures shows the undistorted lattice of N sites with random 
spring constants kg and constant lattice spacing a. The bot- 
tom figures shows that stretched lattice with length -yNa and 
random lattice spacings ag — 7 + u' t — u' l _ 1 . 



where a is the rest bond length. Harmonic springs with 
spring constant kg = k + 5kg connect sites £ and I — 1, 
where k = (J2 g kg)/N is the average spring constant and 
J2e 5kg — 0. The lattice is stretched from its equilib- 
rium length Na to a new length L = jNa. If all kg's 
were equal, the lattice would undergo an affine distor- 
tion with Rg = ja£. When the k/s are random, sites 
undergo an additional nonaffine displacement u' e so that 
R e = jai + u'f. The energy is thus 



1 N 



(3.6) 



In equilibrium, the force Fg — ~dH./du' e on each bond is 
zero. The resulting equation for u' £ is 

Fg = ke + i('ya+u'£ +1 -u' i )-ke('ya+u't-u't_ l ) = 0, (3.7) 

which can be rewritten as 

-A+fc £ A_z4 = jaA + Sk h (3.8) 

where A + and A_ are difference operators defined via 
A+Ag = Ag + \ — Ag and A_ = Ag — Ag-i for any function 
Ag. The Fourier transforms of A + and A_ are, respec- 
tively, A+(q) = e lq — 1 and A_(g) = 1 — e~ lq . Equations 
(|3.7(l and 1)3. 8|1 must be supplemented with boundary con- 
ditions. We use periodic boundary conditions for which 



u' or equivalently 



N 

1=1 



(3.9) 



The solution to Eq. 1)3. 8|1 can be written as the sum of a 
solution, 



-(A + foA_) 1 jaA + 5kg = -jaA 



_j 5k e 



(3.10) 



to the inhomogeneous equation and a solution, 

A+k e A-uf = 0, (3.11) 

to the homogeneous one. The latter solution is uf = 
A_C/kz where C is an as yet undetermined constant. 
Adding the two solutions we obtain 

1 (3.12) 



^ a ki + kg 



which implies A-u' e — (~ja5kg + C)/kg. The boundary 
condition of Eq. 13. 9|) determines C, and the final solution 
for u'g is 

u' t = - 7 aA: 1 ^- skt - (J2 K 1 )^ 1 K l5k t 



= 7aA_ 1 



— -11= -jaAZ 1 S i . (3.13) 



The quantity 



S. = l- 



{i+nY 



(3.14) 



depends only on the ratio pg = 5kg/ k. 

Equation (|3.13|) is the complete solution for u'g for an 
arbitrary set of spring constants kg. In our model, these 
spring constants are taken to be random variables, and 
information about the nonaffinity is best represented by 
correlation functions of the nonaffine displacement, av- 
eraged over the ensemble of random kg's. The simplest 
of these is the two-point function G{£ — £') — (ug,u' e ), 
where ( } represents an average over kg. G(£ — £') is eas- 
ily calculated from Eq. 1)3.13)1 : its Fourier transform is 



G(q) = ( 7 a) 2 



A s (q) 



(3.15) 



2(1 — cosq) 

where A s (q) is the Fourier transform of A s (£' — £) = 
(Sg'Sg). 

There are several important observations that follow 
from the expression Eq. I|3.15|l and that generalize to 
higher dimensions. 

• A s (£',£) depends only on the ratios 5kg/k and 
Skgi /k, and it increases with increasing width of the 
distribution of 5kg. To lowest order in averages in 5kg, 

A s (£\£) = k- 2 [A k {£',£)-N- iy £ j A k {i,£ 1 )} 



k~ 2 A k (5g,g, -N- 1 ), 



(3.16) 



where A k (£',£) = (5kg' 5kg). The final form applies to mi- 
correlated distributions in which spring constants on dif- 
ferent bonds are independent and (SkgSkgi) = A k 5g.gi . As 
the width of the distribution increases, higher moments 
in 5kg become important in A s (£',£). If we assume that 
the only nonvanishing fourth order moments are of the 
form A^ (£',£) = ((5k t ) 2 (5kg) 2 ), then the fourth-order 
contributions to A s are 

fc 4 A^ 4 > (£',£) 



= 1- 



4 6 
N + 



) A^(£',£)-^Y, AikA) (^^ 
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For uncorrelated distributions, A^"' 4 ^ (£, £') = 
A^Sej, + (A fc ) 2 (l - 6 e<e ,). Thus, for uncorrelated 
distributions in the limit TV — > oo, 



A k 3A^ 



(A fc ) 2 



fc 4 



(l-(5,, ) (3.18) 



to fourth order in Skg. Note that the constraint Skg = 

requires Y,e A 5 (^, e ) = £f A 5 (^, l ') = and thus that 
A s (q = 0) = 0. This condition is imposed by the factor 

1 — <5 g ,o in Eq. l|3.18[l . which implies that lim g ^o A s (q) ^ 
A s (q = 0) = 0, i.e., A s (q) does not approach zero as 
some power of q as q — > 0. 

It is easy to verify that Eq. l|3.18[l is exactly the same 
result that would have been obtained using only the 
solution [Eq. (|3.10[l ] to the inhomogeneous equation for 
u' e with Ski/ki replaced by 5kg /kg — N~ 1 '^ li Skg/kg. 
Thus, to obtain the solution for G(q) to leading order 
1/N, we can ignore the boundary condition, Eq. I|3.9fl . 
and use the solution to the inhomogeneous equation 
with the constraint that q 2 G(q) be zero at q = 0. This 
observation will considerably simplify our analysis of the 
more complicated higher-dimensional problem. 

• If correlations in Skg are of finite range, then 
A s (q) has a well defined q — > limit. In this limit, 

G(q) = (7«) 2 ^ « {J S^(1 = 0), (3-19) 

where A A (Q) = lim^o A A (q) for A = S,k. Thus, there 
is a q~ 2 divergence in G(q), and the spatial correlation 
function Q(£',£) = ((u' e — ug) 2 ) diverges linearly in sepa- 
ration 

/\ k ((\\ 

Q{t,e) ~ ( ia ) 2 A s (0)\£'-£\ ~ ( ia ) 2 —^\£'~£\. (3.20) 



• If correlations in Sg extend out to a distance £, 
then A s (q£) becomes a function of q£. Long-range 
correlations in Sg will lead to long range correlations in 
G(q) ~ q~ 2 A s (q£), and Q{£, 0) will grow more rapidly 
than £ for 1 <C £ <C £• It is possible that this is the 
correlation length that diverges at the jamming point 
J in granular media [35L |3(j . We will discuss this point 
further in Sec. liTTEl ' 



B. Nonafflnity for d > 1 

The nonaffinity correlation function, 

G ij (x,x') = «(x)n;(x')) ! (3.21) 

for d > 1 has a form very similar to that for d = 1, 
except that it has more complex tensor indices. We will 
be primarily interested in the scalar part of this function, 



obtained by tracing over the indices i and j. The Fourier 
transform of this function scales as 



G(q) = G»(q) ~ ^A S (q) 



7 2 A*(q) 



K 2 



(3.22) 



where 7 represents the appropriate components of the ap- 
plied strain and A (x, x') is in general a nonlinear func- 
tion of the ratio of the fluctuating components SKijki(x) 
of the elastic-modulus tensor to its uniform components 
Kijki- To lowest order in the variance, A s ~ A K / K 2 
where A^ represents components of the variance of the 
elastic-modulus tensor and K components of its average. 
Thus, the nonaffinity correlation function in coordinate 
space is proportional to |x| — ( d ~ 2 ) in dimension d, or 



((u'(x)-u'(O)) 2 } 
Aln(|x|/B) d=2 
C-DM- 1 d=3, 



where 



.4 
B 

C 

D 



- 7 2 A S (0) p 
(aA)- 1 

7 2 A S (0)A 

- 7 2 A S (0), 

7T 



1^A*(0) 

7T 



K 2 



(3.23a) 
(3.23b) 
(3.23c) 



(3.24a) 
(3.24b) 
(3.24c) 

(3.24d) 



where A = 2n/a is the upper momentum cutoff for a 
spherical Brillouin zone with a the short distance cutoff 
and a — 0.8905 is evaluated in App. [U] The length B 
depends on the spatial form and range £ of local elastic- 
modulus correlations. We will derive explicit forms for it 
shortly. In our numerical simulations, we allow the bond 
spring constant k^ to be a random variable with variance 
A fc = ((Ski,) 2 ). Variations in kb in general induce changes 
in all of the components of SK^i , and A is an average 
of a function of Sk^/k where k is the average of fcf,. 

In general <7(x) also has anisotropic contributions 
whose angular average is zero. We will not consider these 
contributions in detail, but we do evaluate them analyt- 
ically in App. IU1 

When a sample is subjected to a distortion via stresses 
at its boundaries, the strains can be expressed in terms 
of an affine strain and deviations from it. Using the ex- 
pressions in Eq. (|3.4() for these strains, we obtain the 
energy 



SH = 



\ J {Kijktdju'idtu'k 

+ [5K ljk i(x.) + Sik&jiix^djUidiu'k 
+25K ijH (^) lM d j u' i } (3.25) 



to lowest order in 7^ . 
we obtain 



Minimizing STL with respect to u', 



-dj [K ijU + SK ijU (x) + S ik aji (x)]cWfc 



djSK 



ijkl 



( x h 

(3.26) 



A:/- 
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This equation shows that the random part of the elastic- 
modulus tensor times the affine strain acts as a source 
that drives nonaffine distortions. The random stress, 
which is transverse, does not drive nonaffinity; it is the 
continuum limit of the random force. The operator 
-djK^WdtSix-x') ee x^x'), where ^„(x) = 
Kijki+3Kijki{it) + Sik<7ji(it) is the continuum limit of the 
dynamical matrix or Hessian discussed in Refs. [25j and 
|24|. The matrix Xij(x, x') is the response Sui(yi)/Sfi(x') 
of the displacement to an external force. The formal so- 
lution to Eq. (|3.26[1 for u\ (x) in terms of SK^i (x) and 
&ij (x) is trivially obtained by operating on both sides 
with x P i( x > x ') : 



uj(x) = J d d x' Xip (*-^)d' J 5K nkl {x!) lku (3.27) 

The random component of the elastic modulus appears 
both explicitly and in a hidden form in \ip m this equa- 
tion. 

Equation l|3.27[) is the solution to the inhomogeneous 
equation, Eq. (|3.26|) . Solutions to the homogeneous equa- 
tion should be added to Eq. (|3.27() to ensure that the 
boundary condition u'(xb) = for points xj on the 
sample boundary is met. As in the ID case, however, the 
contribution from the homogeneous solution vanishes in 
the infinite volume limit and can be ignored. 

To lowest order in the randomness, we replace \ip m 
Eq. (|3.27() with its nonrandom counterpart, Xi P ( x — x ')> 
the harmonic elastic response function <5iii(x)/£/j(x') of 
a spatially uniform system with clastic-modulus tensor 
Kijki to an external force /j(x'). Thus, to lowest order 
in 7«) 

GWfa) = X°p(q)x°p'(-q)9j*'A^ fei;p , /fc , i ,(q)7fci7fc'i', 

(3.28) 

where 



A^fci;i'j'fe'i'( X 7 X ') 



(SKijuWdKi^kn,^')) (3.29) 



is the variance of the elastic-modulus tensor, which we 
simply call the modulus correlator. Equation (|3.28(l 
contains all relevant information about nonaffine corre- 
lations to lowest order in the imposed strain. It applies 
to any system with random elastic moduli and stresses 
regardless of the symmetry of its average macroscopic 
state. 

Our primary interest is in systems whose elastic- 
modulus tensor is macroscopically isotropic. In these 
systems, which include two-dimensional hexagonal lat- 
tices, K^m = XSi 3 d k i + fJ-iSikdji +6u6jk) is characterized 
by only two elastic moduli, the shear modulus /i and the 
bulk modulus B = A+(2/i/d), where d is the dimension of 
the reference space. The Fourier transform of x^ (x, x') 
in an isotropic system is 



4(q) 



i 



i 



(A + 2n)q< 



9iQj H ~ ( 3 - 30 ) 



The modulus correlator is an eighth-rank tensor. At 
q = 0, it has eight independent components in an 



isotropic medium (See App. ^ and more in media with 
lower symmetry, including hexagonal symmetry. As dis- 
cussed above, however, all components of SKijki are pro- 
portional to Skb. 

We show in App. [B]that G(q) has the general form 



G(q) 



7 



■<•</ 



[i 2 q 2 



(A A + A B q 2 ± - Acfjl), 



(3.31) 



where q t = qi/q, q\ = q% + q 2 . and A^, A B and A c 
are linear combinations of the independent components 
of Af- kli , j, kn , times a function of X/^i. Thus, in general 
(?(x) will have anisotropic parts that depend on the di- 
rection of x in addition to an isotropic part that depends 
only on the magnitude of x. In App. [U] we derive ex- 
pressions for the full form of <?(x). Here we discuss only 
the isotropic part, which has the from of Eq. (|3.23() with 



.4 



7z y 
TTfl 2 



[A, 



1 



B = (aAY 



C 



D = 



i X 



7ai/ 

27T/i : 



[A, 



■[A, 



-A, 



15 



A C ]A 



2 A 1 A 1 

3 - 15 Ad 



(3.32) 
(3.33) 
(3.34) 

(3.35) 



In two dimensions, the anisotropic term is proportional 
to cos 4-0 where t/j is the angle that x makes with the 
x-axis. In the limit of large |x|, the coefficient of cosAip 
is a constant. In three dimensions, the anisotropic terms 
are more complicated. In both two and three dimensions, 
however, the average of the anisotropic terms over angles 
are zero. 



C. Other Measures of Nonafnnity 

The nonaffinity correlation function G^ (and its cousin 
Q) \s not the only measures of nonafinity, though other 
measures can usually be represented in terms of it. 
Perhaps the simplest measure of nonaffinity is simply 
the mean-square fluctuation in the local value of of 
u'(x), which is the equal-argument limit of the trace of 
Gij(x',x): 



<[ U '(x)] 2 )=G l4 (x,x). 



(3.36) 



This measure was used in Ref. |8j to measure nonaffin- 
ity in models for foams. In three dimensions, it is a 
number that depends on the cutoff, a : ([u'(x)} 2 ) ~ 
j 2 (A K /K^a^ 1 ; in two dimensions, it diverges logarith- 
mically with the size of the sample L: ([u'(x)] 2 ) ~ 
7 2 (A X /K 2 )\n(L/al 

References |2(j, [^J [2^, which investigate a two- 
dimensional model of crosslinked semi-flexible rods de- 
signed to describe crosslinked networks of actin and other 
semi-flexible biopolymers, introduce [Fig. [3] a measure 
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'(x') X 



x',x 



R^(x') 



R(xL 



u'(x) 



e (x-,x) 



R A (x) 

FIG. 3: Graphical representation of the angular measure of 
nonaffinity used in Ref. |20ll . Under afEne distortions, points 
Xi transform to points i?Ai( x ) = x i+lij x j> an( A the vector r = 
Ra(x') — Ra(x) makes and angle #o(x',x) with the x-axis. 
Under nonaffine distortions, points Xi transform to -Ri(x) = 
_Rai(x) + iti(x), and the vector r' = R(x') — R(x) makes an 
angle #(x',x) with the i-axis. 



based on comparing the angle 9 = #(x',x) that the vec- 
tor connecting two sites originally at x and x' makes with 
some fixed axis after nonaffine distortion under shear to 
the angle 8q = 6>o(x') — Oq(x) that that vector would make 
if the points were affinely distorted: 



<5 e (x'-x) = <[0(x\x)-# o (x',x)] 2 ). 



(3.37) 



Under affine distortion, the vector connecting points x' 
and x is n — x[ — Xi + ^ijix'j — Xj); under nonaffine 
distortion, the separation is r' = r + u'(x') — u'(x). In 
two dimensions, 



r x r = rr sin 



o )e z = rx [u'(x') - u'(x)], (3.38) 



where e z is the unit vector along the z direction perpen- 
dicular to the two-dimensional plane and r = |r|. If both 
7 and |u'(x') — u'(x)|/|x' — x| are small, 



0(x',x)-fl o (x',x) 



and 



g (x) = eye fe i^f£fcz(x) 



e z -[(x'-x)x(u'(x')-u'(x))] 



(3.39) 



(3.40) 



where = e Z ij is the two-dimensional antisymmetric 
symbol, and ft,(x) = <[u{(x) - uj(0)][«^(x) - u$(0)]). 

D. Generation of Random Stresses 

As we have discussed, a system of particles in mechan- 
ical equilibrium can be characterized by random elas- 
tic moduli and a random local stress tensor with only 



transverse components. To better understand random 
stresses, it is useful to consider a model in which ran- 
dom stress is introduced in a material that is initially 
stress free. We begin with a system with a local elastic- 
modulus tensor Kijki(x) that can in general be random 
but with <7jj (x) = 0, and to this we add a local random 
stress Uij (x) with zero mean that couples to the rotation- 
ally invariant nonlinear strain and that has longitudinal 
components so that its variance in an isotropic system is 



(3.41) 



For simplicity, we assume that the spatial average of 
Wij (x) is zero. A random stress of this sort can be gen- 
erated in a lattice model by making the rest bond length 
RbR a random variable in a system in which initially the 
rest and equilibrium bond lengths are equal. In the con- 
tinuum limit, our elastic energy is thus 



8H 



[5^h( x KW°h( x )+^( x K( x )]. ( 3 - 42 ) 



where the superscipt a on K^ kl indicates that this is an 
clastic modulus prior to relaxation in the presence of <zy . 

Sites that were in equilibrium at positions x in the 
original reference space in the absence of TTij are no longer 
so in its presence. These sites will undergo displacements 
to new equilibrium sites x' = Ro(x) = x + Uo(x), which 
define a new reference space. Positions R(x) in the target 
space can be expressed as displacements relative to the 
new reference space: R(x') = x' + u'(x'). Then, strains 
relative to the original reference space can be expressed as 
the sum of a strain relative to the new reference space and 
one describing the distortion of the original references 
space to the new one: 



Uy(x) 



1 fdR k (x) dR k (x) 



2 \ dxi dxj v 



where 



u2,(x) + A^x'y^xOAoytx'), (3.43) 



A y(x) = ^*W =Jy+ iUoii (3.44) 



dXj 

U i] = 2 (^O/ciAofcj - Sij) , 



(3.45) 



and 



where d[ = d/dx^. Using Eq. lpH3l) in Eq. % J >A'l). we 
obtain SH[u] = 5H[u } + 6H'[u'], where 



< (x') = + dX + dlu'M^/2, (3.46) 



5H' = ~ I ^1%^ + ajiu'ij] (3.47) 



with 



K m (x') = (det A r 1 Ao ia Ao ]b K bcd Al k ^odi (3-48) 
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and 



E. Long-range Correlations in Elastic Moduli 



c^-(x') = (detA ) 1 A 0t a(KbcdUocd + Vab)Aob J , (3-49) 

where we have not displayed explicitly the dependence of 
Aoij on x' . The displacement field Uo (x) is determined by 
the condition that the force density at each point in the 
new reference state be zero, i.e., so that d'j&ij(~x!) = 0. 
To linear order in displacement and Wij , this condition is 



dj(K ijk iu ok i + (Tij) = 0, 



(3.50) 



where to this linearized order, we can ignore the differ- 
ence between x and x'. For an initially isotropic medium, 
this equation can be solved for Uo to yield 

i \ 1 fx A + V 9iQk\ ■ - fQK1 s 
fxq 2 v A + q 2 j 

To lowest order in uo, the elastic moduli and stress ten- 
sors in the new reference state are 



^y-(q) = 8iktfi<rki - x + 2it 5 ' i ™ l ' ru ( ' ! ' r>21 



SKijM = 2X(5ijVki + 6kmj) (3.53) 
+2fi(6ikVji + Sjmk + SuVjk + SjkVu), 



where 



2 (AoifcAojy - hj) 

(diu 0j + djU 0i + d k u 0i d k u 0j )/2 

{diu 0j + djU i)/2 



(3.54) 



is the left Cauchy strain tensor relative to the original 
reference state. 

Note that ov,(x') is transverse and random as it should 
be. The elastic-modulus tensor is a random variable via 
its dependence on Aoiy(x). Thus, a random stress added 
to an initially homogeneous elastic medium (with Kfj kl 
nonrandom and independent of x) produces both a ran- 
dom transverse stress and a random elastic-modulus ten- 
sor in the new relaxed reference frame. The statistical 
properties of Kij k i(x.') are determined in this model en- 
tirely by those of <7y(x), and A K ~ A a . In general, 
of course, the randomness in Kij k i(x') arises both from 
randomness in the original Kfj kl (pc) and ffij(x). 

The nonaffinity correlation function can be calculated 
exactly to lowest order in and A£ when the initial 
reference state is homogeneous and nonrandom. It has 
exactly the same form as Eq. (|3.22|l when expressed in 
terms of A^. When expressed in terms of A \ and Af , it 
has a similar form, which in an isotopic elastic medium 
can be expressed as 

G(q) ~ ^/(q, A//i, A?/Af). (3.55) 
Thus, <?(x) has the same form in this model as Eq. I|3.23[l . 



Long-range correlations in random elastic moduli can 
significantly modify the behavior of <?(x). To illustrate 
this, we consider a simple scaling form for A K (q) inspired 
by critical phenomena: 



A A (q)-eW) 
^ l><?o[l + (<?0 S + •••], 



(3.56) 
(3.57) 



where £ is a correlation length, <f> is the dominant crit- 
ical exponent, and s and t are corrections to scaling 
exponents. It is possible in principle for each of the 
components of Af- kli , - lk , v to be described by difference 
scaling lengths £ and functions g(u). We will assume, 
however, that £ and the functional form of g is the 
same for all components, but we will allow for the zero- 
argument value go to vary. G(q) is thus given by Eq. 
(|3.31|) with A a , A B , and Ac replaced by A^ (q) , A B (q) , 
and Ac(q) with scaling forms given by Eq. I|3.56[l . In 
this case, £7(x) can be written as (7^ y //x 2 )J r (x) with 
.F(x) = ^a(x) + .Fb(x) - ^c(x), where 



d d q 1 

7^p/a(q)ffa(^)^(l 



(3.58) 



with f A = 1, /s(q) = eg. , and / c (q) = g 2 ^ 2 ,. There are 
two important observations to make about the functions 
T a . First, for q£ <C 1, g(q£) can be replaced by its zero q 
limit, g . Thus, as long as £ is not infinite, the asymptotic 
behavior of Q(x) for |x| 3> £ is identical to those of Eq. 
(|3.23|l but with amplitudes that increase as Second, 
when £ — > oo, the g( d_3 ~^) behavior of the integrand 
leads to modified power-law behavior in |x| for a <C x -C 
£, where a = 2tt/A is the short distance cutoff, depending 
on dimension. 

In two dimensions, which is the focus of most of our 
simulations, the isotropic part of T is 

^(x) = £. [ %(gO[l - J (ff|x|)], (3.59) 

where g(y) = g A {y) + gB{y) - \gc{y) and J (y) is the 
zeroth order Bessel function. In the limit |x| 3> £, 



i KH^) 

J>(x) 5oC^ln 7 x, 

7T £ 



(3.60) 



where /3(A£, cj>) is evaluated in Appendix[U] The behavior 
of .Fj(x) when A -1 <C |x| <C £ depends on the value of (f> 



if 6 < 2, 



^(x)~<(^ ffco |x| 2 lnK/|x|), if = 2, 
^5o O |x| 2 ^- 2 C 2 (0), if > 2. 



(3.61) 



The quantities A2(4')i ^(0), and ^ are evaluated in Ap- 
pendix |2| 
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FIG. 4: The function (5(|x|) in two dimensions for fixed 7 and 
K for systems with long-range correlations in A K defined by 
Eqs. 13.591 and 13. 631 for different exponents (f>. The ampli- 
tudes and correlation lengths £ for each <f> are normalized so 
that all curves have the same values of B = £//3(A£,0) and 
coefficients of In |x|/_B at large |x|/B. The curve for (f> — 2 
corresponds to Lorentzian correlations [Eq. 13.621 1. [Color 
online] 



The function g{u) can have any form provided its large- 
and small- it limits are given by Eq. (|3.56|) . A useful model 
form to consider, of course, is the simple Lorentzian for 
which d> = 2 and 



(3.62) 



For the purposes of illustration, in Fig. 0]we plot JFj(|x[) 
for a family of functions parameterized by the exponent 



9(u) = 



9o (</>) 



(l + u 2)0/2- 



(3.63) 



these curves clearly show the crossover from |x|* behavior 
for A -1 <C |x| <C £ to the characteristic log behavior for 
x| » £. The correlation length £ and the amplitude 
go((j)) were set so that the large x log behavior is the 
same for every For this family of crossover functions, 
the value of |x| at which |x| ) crosses over from |x|^ 
to logarithmic behavior increases with decreasing 0, and 
curves with smaller (f> systematically lie above those with 
large cf). 

The limiting forms for .F/(|x|) in one and three dimen- 
sions are given in App. [U] 



F. Rotational Correlations 

The nonaffinc displacements generated in random elas- 
tic media by external strains contain rotational as well 
as irrotational components as is evident from Fig. JSJ. 
The local nonaffine rotation angle is Wfe(x) = ^€ijkdju' k , 
where e^fe is the anti-symmetric Levi-Civita tensor, and 



rotational correlations are measured by the correlation 
function G u , iUJj (x) = (ti>;(x)cjj(0)}. In two dimensions, 
there is only one angle uj(x) — ^e r id r Ui, where e r i = 
e zr i. The Fourier transform of the correlation function 
G u = (w(x)w(O)) will then scale as 7 2 A K //i 2 , approach- 
ing a constant rather than diverging as q — * 0. We show 
in App. |D] that 



GUq) 



7 2 
M 2 



[A» A (q) - A» c (q)q 2 J 



(3.64) 



are linear com- 

K 



in two dimensions, where A^(g) and A£ 
binations of the independent components of A.^^, ■, k , t , 
Thus, the rotation correlation function contains direct 
information about elastic-modulus correlations. If these 
correlations are short range, and there is no q depen- 
dence in either A^(g) or A£>, the spatial correlation func- 
tion has an isotropic short-range part and an anisotropic 
power-law part: 



CI 



, x = 



A^(x) - A^ 



16 



2 2 
x y 



(3.65) 

If there are long-range correlations in the elastic moduli 
with the Lorentzian form of Eq. (|3.62l) , then 



G„(x) 



Ixy 



--cos4V>Ag 



*o(|x|/0 



(3.66) 



48f 



4£2 



where K n (y) is the Bessel function of imaginary ar- 
gument. The cos 4^ behavior is for isotropic systems. 
There will be cos6"0 and higher order terms present in 
a hexagonal lattice. In Sec. IIVI we verify in numerical 
simulations the exponential decay of the isotropic part 
of G w (|x|) in Model A with long-range correlations in 
spring constants and the |x| -2 behavior of the cos 4-0 
part of G w (|x|) in Model C, which is isotropic. 



IV. NUMERICAL MINIMIZATIONS 

To further our understanding of nonaffinity in random 
lattices and to verify our analytic predictions about them, 
we carried out a series of numerical studies on models A- 
D described in Sec. Ill Al To carry out these studies, we 
began with an initial lattice - a periodic hexagonal or 
FCC lattice for models A and B and a randomly tesal- 
lated lattice for models C and D. We assigned spring po- 
tentials Vb(Rb) and rest bond lengths RbR to each bond. 
To study nonaffinity, we subjected lattices to shear and 
then numerically determined the minimum-energy posi- 
tions of all sites subject to periodic boundary condition. 
The elastic energy of the lattice was linearized about the 
affine shear state. Interestingly, in this linearization the 
value of the imposed shear, 7, factored out of our cal- 
culation, so that u'(x) was linear in 7 and thus S(x) 
was automatically quadratic in 7. We present below the 
procedures and results for each model. 
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FIG. 5: This figure shows the direction of nonaffine displace- 
ments resulting from a shear in the xy-p\ane displayed as 
arrows on the original reference lattice. Note the vortex-like 
patterns. 



A. Model A 

In this model, the initial reference lattice is periodic, 
and the rest bond length R b R is equal to the equilibrium 
lattice parameter R b Q for every bond, which we set equal 
to one. Each bond is assigned an anharmonic potential 

V b (R b ) = h b (8R 2 + 8R 4 b ), (4.1) 

where 8Rr = R b — R b Q = R b — 1 and the spring constant 
k b is a random variable. We chose k b = 1 + 6k b where 
5k b is a random variable with zero mean lying between 
— 8k and +8k with 8k < 1. 



1. Independent bonds on hexagonal and FCC lattice 

In the simplest versions of model A, the spring constant 
k b is an independent random variable on each bond of a 
two-dimensional hexagonal or a three-dimensional FCC 
lattice. We assign each bond a random value of 8k b cho- 
sen from a flat distribution lying between — 8k and +8k. 
Randomly distributed spring constants give rise to ran- 
dom local elastic moduli as defined by Eq. H2.17fl . We 
verified that the distribution of the values of the local 
shear modulus K xyxy on a hexagonal lattice for different 
8k was well fit by a Gaussian function with width linearly 
proportional to 8k. 

The nonaffinity correlation function £?(x) = 
(|u'(x) - u'(0)| 2 ) [Eq. lEH] measured on the nu- 
merically relaxed lattices is shown in Fig. HJa). The 
averages were calculated by summing the differences 
in deviation for every pair of nodes on the lattice 
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FIG. 6: The nonaffinity correlation function Q{x)j ('ySk) 2 for 
sheared model-A lattice: (a) Q(x)/(j8k) 2 versus |x| for a 2D 
hexagonal lattice sheared to 0.1% for values 8k = O.Oln for 
n = 1, 10 of the variation in local spring constant, (b) 
g{x)/{8k) 2 versus |x| for a 3D FCC lattice sheared to 0.1% 
for several different values 5k = O.Oln for n = 1, 10. The 
data is fit to a function C — _D/|x|. All lengths are in units of 
the lattice spacing. [Color online] 



and binning according to the nodes' separation in the 
undeformed (reference) state. Note that this process 
automatically averages over angle, so it produces only 
the isotropic part of £/(x). The separation between 
nodes was taken as the least distance between the nodes 
across any periodic boundaries. The curves were well 
fit by the A\n(\x\/B) dependence on |x| predicted by 
Eq. (|3.23b|) . The excellent data collapse achieved by 
plotting the rescaled function <5(x) / (jSk) 2 demonstrates 
the quadratic dependence of the amplitude A on 8k. 
Figure [7| shows the quadratic plus quartic dependence 
of the amplitude A on and 8k at larger values. It 
is worth noting that while all correlation functions 
were independently fit with a two-parameter function 
Aln(|x|/£>), the optimal values of B in all cases fell 
within 10% of one another. 

FigureHJb) displays Q(x) / (8k) 2 on an FCC lattice as 
a function of |x| for different 8k for 7 = 0.1% fit to the 
function C - D/\x\ predicted by Eq. (|3.23|l . The data 
collapse verifies the expected dependence of C on (8k) 2 . 
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FIG. 7: (a)The prefactor A of the fits in plots of Fig. EJa) 
as a function of <5fc. The dotted curve is a fit to a quartic 



g(Sk) 2 + h(8k) 4 as suggested by Eq. (JHUHl ■ 
data up to 5k = 0.1 and a quadratic fit. 



The inset shows 



o 




FIG. 8: The two point correlation function A K (x = 0) = 
{6K X y X y(x)6Kxyxy(0)) produced by random spring constants 
drawn from an exponentially correlated random scalar field. 
The exponential decay is evident. [Color online] 



2. Correlated random bonds on an hexagonal lattice 

As discussed in Sec. IIII El random lattices can ex- 
hibit long-range correlations, characterized by a corre- 
lation length £, in local elastic moduli that can signif- 
icantly modify the behavior of nonaffinity correlation 
functions at distances less than £. To verify the predic- 
tion of Sec. IIII Dl we numerically constructed hexagonal 
lattices with long-range correlations in bond spring con- 
stants. To do this, we set kf, = 1 + Skb where 5kb was 
set equal to a small, randomly generated scalar field with 
proper spatial correlations. This scalar field was created 
by taking the reverse Fourier transform of the function 
exp(i(f> r ) I \J q 2 + £~ 2 , where £ is a variable decay length 
and 4>r is a random complex phase. The scalar field in 
these simulations was normalized to have constant mean 
squared value and peak values of ±0.1, so that the vari- 
ation to the local spring constants was at most 10%. 
This method of generation yields a clean exponential 
decay in the two-point correlation function A (x, 0) = 
(SK xyxy (x)SK xyxy (0)) which persists for separations up 
to three times the correlation length. Figure |H1 shows the 
two-point correlation function A K (x, 0) as a function of 
separation. The region of exponential correlation was fol- 
lowed by a small region of anti-correlation, which is not 
pictured. By construction, the distributions of the local 
shear elastic modulus K X y Xy were essentially constant, 
independent of £; thus A (0, 0) is equal for all curves in 

Fig. m 

According to Eq. (|C13I) . the growth of the correlation 
function <5(x) for large |x| is logarithmic with prefactor 
proportional to <7o£ j where = 2 for the Lorentzian case 
we are now considering. The quantity go^ is equivalent 
to A x (q = 0), but this quantity is difficult to measure 
numerically. However, 170 can also be expressed in terms 
of the coordinate space correlation A K (x = 0). The latter 

quantity is easily measured by averaging (^(5K xyxy (x)) 2 ^ 

over all nodes. For the form of the correlation function 



g{u) given in Eq. i|3.63[l . 

^=o,~. w ^H i+ < A{,2 r*) 

[iln(l + (A0 2 ) = 2. 

(4.2) 

In the limit £A 3> 1, go ~ (x = 0) and the large sepa- 
ration form of the correlation function Q (x) is logarithmic 
with prefactor 7 2 A K (x = 0)^//i 2 . 

We have already established that for this set of sim- 
ulations, A K (x = 0)/fi 2 is a constant, independent of £ 
[see Fig E| . In Fig. E3 we plot £7(x)/£ 2 versus |x|/£ for 
different values of £. We also plot the function .F(|x|) 
calculated from Eq. (|3.59l) with a Lorentzian g(y) [Eq. 
(|3.62l) ]. The agreement between the numerical and ana- 
lytical results is excellent with both showing |x| behavior 
for |x| < £ and ln|x| behavior for |x| > £. In Fig. 1101 we 
plot the vorticity correlation function G u versus separa- 
tion rescaled by the correlation length, |x|/£. The vor- 
ticity correlation function decreases exponentially away 
from zero separation with a decay length ~ 1.1 x £; our 
framework predicted decay with an exponent of £ exactly. 
The slight discrepency between theory and simulation is 
not understood. 



B. Model B: Internal stresses 

In this model, random stresses are introduced in a 
periodic lattice via a random distribution of rest bond 
lengths. We study hexagonal lattices in which the rest 
lengths of the bonds are multiplied by a factor (1 + fib) 
where (3b is chosen randomly from a flat distribution ly- 
ing between —(3 and (3 with (3 < 0.1. Once again, the 
spring constants are set to kb = 1 + 5kb, with Skb cho- 
sen randomly from a flat distribution lying between — 5k 
and +5k. After specifying the rest length of each bond, 
we numerically determined the equilibrium state of this 
random lattice with zero applied stress by minimizing 
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FIG. 9: Plot of G(x)/£, 2 versus |x|/£ for different £ from 
(points) numerical minimizations and (solid line) the analyt- 
ical expression of Eq. HM.59|I with a Lorentzian g(y). [Color 
online] 
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FIG. 11: The variance of the local shear modulus versus j3 for 
several different 8k. [Color online] 
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FIG. 12: J 4/(7 2 A K (x = 0)/^i 2 ) versus f3 for several different 
Sk. [Color online] 



FIG. 10: Plot of the vorticity correlation versus |x|/£ for 
several different £ for lattices with Lorentzian spatial correla- 
tions in Ski,. [Color online] 

the rest energy over lattice positions and the size of the 
simulation box (for a system of 40,000 particles, the min- 
imization over box size was only a fraction of a percent). 
The resulting equilibrium configuration has zero net force 
on each node. This relaxed state constitutes the refer- 
ence state of our random system with lattice positions 
R<?o = x. 

The original lattice before relaxation is characterized 
by random stresses cFy , which can be calculated from Eq. 

W = Y Yj RbiiRbijk b 5R h /a, (4.3) 
v 

where Hm is the bond vector of length a (independent of 
b) for bond b in the initial undistorted hexagonal lattice 
and SRb = a — RbR = fya. The average over of ay over 
(3b is zero: (&ij)b = 0, and its variance is 

(3 2 (fk) 2 

fiij{£)vki{£))b = -^(1 H —)(SijSki + SikSji + SuSjk). 

(4.4) 

As discussed in Sec. IIII Dl randomness in <7j» generates 
a random elastic moduli in the relaxed reference lattice. 



Figure^Jshows how the random stress broadens the dis- 
tr ibution of l ocal elastic moduli. For lattices with Sk — 0, 
v /A K (x = 0) is linearly proportional to (3 as predicted by 
Eqs. $nn\ to E£5Hl . 

After constructing the relaxed state, we sheared it in 
the xy plane as before and measured the nonaffinity cor- 
relation function. The measurements were well fit by the 
functional form <5(x) ~ Am(|x|/i?). Figure IT21 shows 
that for Sk = the measured ratio A/ (j 2 A K (x = 0) / /i 2 ) 
is nearly independent of /3, as predicted in Section IlII Dl 
For Sk > 0, the ratio A/(j 2 A K (x = O)/^ 2 ) is - 20% 
lower at small f3, but asymptotes to the Sk = value 
as /? increases, approaching the asymptote more quickly 
for smaller Sk. The difference in A/(^f 2 A K (x = 0)//z 2 ) 
between stressed and stress-free lattices is most likely a 
higher order effect due to the breaking of hexagonal sym- 
metry as (3 is increased. 



C. Model C: Random lattice 

In this model, the initial reference lattice is geometri- 
cally random. The rest bond length RbR is equal to the 
equilibrium lattice parameter Rbo for every bond, so that 
the reference state is stress free. Our method of gener- 
ating reference lattices of varying randomness is detailed 
below. Each bond is assigned the anharmonic potential 
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FIG. 13: Section of a random lattice created by triangulating 
a snapshot of a Lennard- Jones gas with T = 8.0 and P = 0.05. 



of Eq. H4.ll) . where the spring constant kb(Rb) = k /Rb 
is a constant per unit length of the rest bond length. 

We use the approach followed in j3lj to generate net- 
works with a tunable degree of randomness. We begin 
by simulating a 2-dimensional gas of 40000 point parti- 
cles interacting through a Lennard- Jones potential. The 
procedure outlined in is used to equilibrate the gas 
at a prescribed temperature and pressure, with periodic 
boundary conditions. The gas is equilibrated for 10000 
time steps, after which the particle configurations are 
sampled every 1000 time steps. In this manner we obtain 
40 uncorrelated configurations of the gas at thirteen dif- 
ferent temperature-pressure combinations, with T = 8.0 
and P = 0.025, 0.05, 0.1, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5 
0.6, 0.7, 0.8, and 1.0, all in units of the Lennard-Jones 
potential. 

We use the particle positions from the snapshots of 
the equilibrated gas as the positions of nodes in our ran- 
dom lattice. Each sampled configuration is rescaled to 
have a box length of 1 on each side. The point configu- 
rations are then tesselated using the Delaunay triangu- 
lation, which places a bond between each node and its 
nearest neighbors. The Delaunay triangulation produces 
networks with an average of 6 bonds per node. A result- 
ing lattice is pictured in Fig. [F3l 

The randomness in local elastic moduli as calculated 
from Eq. (|2.17|l is proportional to the distribution of bond 
lengths and bonds per node. In principle, as we take the 
equilibrium gas pressure to zero, the distribution of bond 
lengths will become completely random. Conversely, as 
we increase the pressure past a critical point the simu- 
lated gas begins to crystalize, forming spatial domains 
of hexagonal order separated by grain boundaries. This 
transition should be marked by a growth in the two-point 
correlation of local shear moduli. We fit the non-affinity 
correlation data for a broad range of pressures which 
cross this transition and compare it to the framework 
developed in previous sections. We used Eq. (|2.17|l to 
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FIG. 14: <5(x) vs |x| for lattices with varying degrees of ge- 
ometrical randomness, scaled to have the same asymptotic 
form for large |x|. The solid line is a fit Eq. 1)3.59^ using g(u) 
of Eq. fifity with cj> = 0.4. [Color online] 




FIG. 15: A/(7 2 A K (x = 0)/^ 2 ) versus B for geometrically 
random lattices, where A and B are the fit parameters to 
the functional form !5(x) = ^41n (|x|/B). The line is a fit to 
these data points using Eq. I|4.5|l . assuming that there are 
long range correlations in the elastic moduli. [Color online] 



calculate A 7< (x = 0)//i 2 for each ensemble of random lat- 
tices, while the crystalline correlation length is fit as an 
unknown. 

The lattice is sheared by 0.1% and the energy is mini- 
mized as a function of node position as before. Figure PHI 
shows the displacement correlation function C?(x) as a 
function of separation |x| for lattices with different de- 
grees of randomness. This correlation function shows the 
same logarithmic growth at large |x| as it does in the ran- 
dom spring constant lattices from the last section. 

We fit the measurements of (?(x) to the functional 
form Aln(|x|/i3) at large |x|. This data is shown in 
Figure ED For the very random lattices generated at 
low Lennard-Jones pressure (T — 8.0, P < 0.3) the val- 
ues of A/(~f 2 A K (x = 0)//x 2 ) and B are nearly constant, 
as our framework predicts for the simple case of delta- 
function spatial correlations. However, lattices created 
at higher pressure values (T — 8.0, P > 0.3) showed 
significant growth of both A/('j 2 A K (-x — Q)/^ 2 ) and B 
with increasing pressure, reaching a saturation point at 
around P — 8.0. Visual inspection of the lattices in ques- 
tion revealed subdomains of hexagonal crystalline order- 
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ing. Long range correlations in the connectivity implies 
long-range correlations in the elastic moduli, so we must 
apply the framework developed in Sec. IIII El and Add. IC ll 
in order to fit the data for partially crystalline lattices. 
Once again, we try the functional form in Eq. (|3.63|) for 
the spatial correlations in the elastic modulus. The nu- 
merical value of the factor go(<j>) can be calculated from 
the measured modulus autocorrelation (x = 0) using 
Eq. @2J|. 

The fitting line in Fig. ^| represents a best fit of both 
the correlation exponent <fi and the cutoff length A -1 to 
the form 



as sugg ested by Eq. (|CT3)> . The best fit was achieved for 
<f> = 0.4 and a cutoff length of A -1 w 1.25 lattice spacings. 
The corresponding analytic form of Q (x) calculated from 
Eq. H3.59fl using g(u) from Eq. I|3.63[) is shown by the 
solid line in Fig. El 

To test the predicted [Eq. I|3.66|l ] cos 4ip anisotropy in 
vorticity correlations, we measured G u (x) as a function 
of the angle x makes with the a;-axis. Figure ITfil shows a 
polar plot of G w (x), which clearly shows cos 4-0 behav- 
ior, and the dependence of the cos 40 term on |x|, which 
shows the expected |x| -2 behavior. 

D. Model D: Random lattice with internal stresses 



(a) 




FIG. 16: (a) Polar plot showing the cos 40 modulation of 
the vorticity correlation function G^x) for |x| = 2, 5, and 
10 lattice spacings. (b) A log-log plot of the coefficient of 
cos 40 in G u (x) showing the expected |x| -2 fall off at large 
x|. [Color online] 



Finally, we simulate the most general model for ran- 
dom lattices, in which the rest bond length R^r is not 
equal to the equilibrium lattice parameter i?f,o, and the 
lattice parameters Rbo along with the number of bonds 
per node are random to within some finite distribution. 
Each bond is assigned the anharmonic potential of Eq. 
(|4.1|) . where the spring constant kb(Rb) = ko/Rb is a con- 
stant per unit length of the rest bond length. We used the 
same geometrically random lattices from Section liV CI as 
staring points, then we add bond length frustration using 
the technique from Section IIVBI We multiply the rest 
lengths of all bonds by a factor (1 + (3b) where (3b is cho- 
sen randomly from a flat distribution lying between —(3 
and (3 with (3 < 0.1. We find the equilibrium configura- 
tion of the lattice by minimizing the elastic energy over 
node positions and box size. We then shear the lattice 
by 0.1%, minimize the energy over node positions, and 
measure the non-affinity correlation function C?(x). 

In all these simulations, the correlation function {?(x) 
was well fit by the functional form A\n (|x[/i3). FigurclT7l 
shows a plot of A/('~f 2 A K (x. = 0)/ii 2 ) for all data sets as 
a function of (3. The data points for (3 — correspond to 
the data from Section HV CI their deviation from the ex- 
pected constancy of Aj ( r y 2 A K (x = 0)/ /z 2 ) was explained 
in that section by the growth of a correlation length scale 




FIG. 17: A/(7 2 A K (x = 0) / u 2 ) versus (3 for geometrically ran- 
dom lattices with internal stresses, where A is the coefficient 
of a logrithmic fit to the non-affinity correlation function 5(x) 
for each lattice and (3 is the bond frustration factor used to in- 
duce internal stresses. The lines are for lattices with eight dif- 
ferent temperature-pressure combinations, with T = 8.0 and, 
from bottom to top P = 0.025, 0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 
and 1.0, all in units of the Lennard- Jones potential. [Color 
online] 
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as the system acquires partial hexagonal crystalline or- 
dering. Here we see that as (5 is increased, the long length 
scale ordering is disrupted by the additional randomness, 
and the ratio Aj (-f 2 A K (x=0) / fJ 2 ) decreases toward the 
value for completely disordered lattices. 

V. SUMMARY AND CONCLUSIONS 

Nonafhne distortions are always present in random 
elastic networks subjected to external stress. In this pa- 
per, using both analytical and numerical techniques, we 
study properties of nonafhnity in these systems mani- 
fested in correlation functions of the deviation, u'(x), 
of local displacements from their affine form. We in- 
troduce four models of random clastic networks with 
random local elastic moduli and possibly local random 
stress arising either from randomness in the form of the 
central force potentials between nearest neighbor sites 
or from random connectivity of the the network. In 
all cases, we show analytically and verify with numer- 
ical simulations that random elastic modulus times im- 
posed strain and not random stress act as sources for non- 
afhne distortions. We calculate the nonamnity displace- 
ment correlation function, <7(x) = ([u'(x) — u'(0)] 2 ), and 
the vorticity correlation function, G u (x) = (w(x)w(O)) 
analytically and verify their form in numerical simula- 
tions for systems with both short- and long-range cor- 
relations in local elastic moduli. We show in particular 
that £(x) - 7 2 (((<5iT) 2 }/ J ft: 2 )ln|x| at large x in two di- 
mensions, where 7 is the imposed strain, K is the average 
of elastic modulus, and ((8K) 2 ) is it variance. 

The formalism we develop is general and should be 
applicable to any elastic system that has a well defined 
average shear modulus. It should provide a basis for 
studying nonafhnity in granular media, foams, networks 
of semi-flexible polymers, and related systems. It should, 
in particular, provide a method of calculating correla- 
tion lengths near percolation-like thresholds such as the 
J-point in jammed systems or the rigidity percolation 
point. We have begun [2i| to use these techniques to 
calculate correlation lengths in the former systems which 
we will eventually compared with those calculated from 
the density of states 1351 [36J and to study nonafhnity in 
networks of semi-flexible polymers [3(j ■ 
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APPENDIX A: PROPERTIES OF THE 
MODULUS CORRELATOR 



The modulus correlator ^^kl-,i'j'k'l'( c * s an ^ n rank 
tensor. The number of its independent components de- 
pends on the symmetry of the reference space. In this 
appendix, we will determine the number and form of 
its independent components at q = (strictly speaking 
q — > 0), or, equivalently, at all q when correlations are 
short range and it is independent of q, when the refer- 
ence space is isotropic. In this case, the general form of 

&W;i>fkj>l< = A ^fei;i'i'fc'i'(q = °) must be constructed 
from products of Kroneker J's that distinctly pair all in- 
dices while respecting all symmetries. 

It is useful to recall how this process is carried out 
for the simpler case of the 4th-rank elastic-modulus ten- 
sor Kijki, which is symmetric under interchange of i and 
j, of k and /, and of the pairs ij and kl. Since any 
index can be paired with any of the remaining three 
and there is only one way to pair the remaining two, 
there are three distinct Kroneker-5 pairings, which we 
will call contractions, of the four indices: SijSki, 5n-5ji, 
and SiiSjk- The first of these satisfies all of the symme- 
try constraints, but the second two do not; their sum, 
however, does. The elastic-modulus tensor, therefore, 
has two independent components in an isotropic medium: 
K^ki = SijSki + ^{S ik Sji + SuSjk)- 

^■fjki-i'j'k'i' 1S symmetric under interchange of i and j, 
k and I, i' and j' , and kl and under the interchange of 
the pairs ij and kl and of the pairs i'j' and k'V; and under 
the interchange of the four-plets ijkl and i'j'k'l'. The 
total number of possible contractions of these 8 indices 
is Nt = 7x5x3x1 = 1 05 because any index can be 
contracted with any of the seven remaining indices, any 
one of the six remaining indices can then be contracted 
with any of the other five remaining, etc. Most of the 
individual realizations of these 105 possible contractions 
will not satisfy symmetry constraints; it is necessary to 
find the linear combinations of them that do. Figure ITSI 
provides a graphical representation of the eight distinct 
contraction groups the sum over whose elements satisfy 
all constraints. The elastic- modulus correlation function 
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where q p = q p /q and 5l ', = 5 PP > - q p q p < . Then 
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with 
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(B2) 



(B3) 



(B4) 



where ^^kl;i'fk'V ^ s defined in Eq. (|A1|) . It is straightfor- 
ward but tedious to calculate S£ and S% from Eq. HAlll . 
The results are 



FIG. 18: This figure provides a graphical representation of 
the groups of contractions that are invariant under the sym- 
metry operations that leave ^^ki;i'j'k'V unchanged. It shows 
8 graphs representing contractions, none of which are trans- 
formed into any of the others under any symmetry operation. 
The symmetry operations applied to each contraction will, 
however, generate other contractions. The number of con- 
tractions for each graph produced by performing all symme- 
try operations on it is indicated below each graph along with 
a representation of the graph in terms of Kronecker S's. 



sf = s? = o 



°3 



S£ = 



AT(qi-me y ), 



St = Si = 0; 
Si = AifJl; 
St = A£4<?i 

Si = AT(2 + dq\ - 16g 2 ), Si = Ai(4ql + 16f x f y ); 
Sf = Af [(d - 1) + q\ AfJ 2 y ]; Si = Af (2 + 4^) 
Sl=Ai{2 + (d-2)qi], Si = Aim. 

(B5) 



APPENDIX C: EVALUATION OF 0(x) 



in an isotropic medium can thus be written as 



ijkl;v j' k'V 



V A Ka 

/ ^ ijkl\i' j'k'l' 



= AiSijSkiSi'j'Sk'i' 

+A 2 6ij8ki(5i>k'Sfi> + Si>i'8j>k>) + prime <-> unprime 
+A 3 (Sik5ji + Sii6jk)(6i>k'8fi> + SiH'Sj'i') 
+A i 8 l j5 i >f5kk'5w + 7 perm. 
+A 5 5ij5i>k>SkfSii' +31 perm. 
+A 6 %<W%'£«' +31 perm. 
+A> r 5 ii ,5 jj ,5kw5iv + 7 perm. 
+A s 5u>6jk>Skji8w + 15 perm. 



where Af^ l . l , fk , v 



AiSijSkiSi'j'Sk'i', etc. 



(Al) 
The first 



three terms in A^ kl . v y k , v describe correlations in the 
isotropic Lame coefficients: Ai = ((<5A) 2 }, A2 = (SXSfi), 
and A3 = ((Sfi) 2 ). The other terms represent fluctua- 
tions away from local isotropy. 



In this appendix, we will evaluate the integrals 
[Eq. (15351) 1 



^^/ Q (q)3^)(l-^ X ) (CI) 



in 2D and 3D that make up the function {?(x), where 
/A(q) = 1, /s(q) = q\, and / c (q) = q*q%. 



1. Two dimensions 



In two dimensions, /s(q) = 1 = /A(q), and 



/ c (q) = sin 2 4> q cos 2 <j> q = -(1 - cos40 s ), (C2) 

where <f> q is the angle between q and the x-axis. Using 
the plane- wave decomposition relation 



e l * x = J (<z|x|) + 2^ cosnej„(g|x|), 



(C3) 



APPENDIX B: EVALUATION OF G(q) 

We outline here the calculation of G(q) to lowest order 
in A K in isotropic systems. We use Eq. (|3.28|) for Gy(q) 
and sum over i = j. Using Eq. (|3.30|l . we find 



x° P (q)x° P '(q) 



1 



1 



(A + 2/x) 2 <r 



-qpip' 



(Bl) 



where J n (x) is the nth order Bessel function, <d = (fiq — ip, 
and ip is the angle between x and the x-axis, and the 
orthogonality relation 



1 

2^ 



, ~ \lS nm cosnip n^Q 
> cos nip cos mfc) = < 't , 

Onm n — 



(C4) 
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we find 



where 



£<t> f A da 

F a = F I \ ga ] = ±- fs o (g0[l-Jb(g|x|)] (C5) fUy. ... f°°dy T ,. nilRO 

Jo 9 lna= / — l-Jo(y) -/ —J (y) = -0.1159 

Jo y 7i y 



for a = A,B and 

•Fc = ^ [sc] - cos 4^ [ffc] , (C6) 

where 

F A [g a ] = - [ A -g a (q£)Hq\x\)- (C7) 



Thus, 



where 



.F(x) = Fj [g] + - cos 4^F A [g c ] , (C8) 



<?(<?£) = 9a{q0 + g B (qS) - pc(Qt). (C9) 

We now evaluate the integrals Fr and Fa in the limits 
|x| > £ > A" 1 and A" 1 < |x| < f. 



(C15) 

and a = 0.8905. When g(y) = goi independent of y, 



ln/?(A£) ^ In a + In AC, 



(C16) 



and Eq. (|C10|) reduces to Eq. H3.23b(l when g^ is iden- 
tified with A s . 

The |x| ^> £ limit of Fa is obtained by setting y = q\x\ 
and noting that letting the upper limit of the integral go 
to infinity and replacing g(y|x|/£) by g Q introduces no 
singularities. The result is 



F A [g]=tg r d Ji M y ) = °£ . 

^ Jo y 4 



(C17) 



a. |x| ^> £ > A 1 in Two Dimensions 



b. A 1 <C |x| <C £ m Two Dimensions 



To evaluate the first limit of Fj, we set y = q\x\ in Eq. 



r /" A|x| dv 
FM = -{ / -3(^/|x|) 

( l ^-g(y^\)[l-J Q (y)]- j * ^-g(y^/\ X \)J (y)} . 

o v Ji y J 



(CIO) 

In the limit |x|/£ — > oo, we can safely replace g(?/£/|x|) 
by g Q in the second and third integrals in this expression, 
and we can let A|x| — » oo in the third integral. The first 
integral diverges as In |x| if we replace g(y£/\x\) by g in 
it, and we have to be more careful to extract the constant 
term beyond the log: 



A|x| 



dy 



du 



g(u) 



«/i> 



At 



du 



g(u) 



f 1 du 

ffo ln|x|/e+ / — [g(u)-g } + 
Jo u 



Thus in the limit |x| » £ > A -1 , 



du 



(Cll) 

g{u)- 

(C12) 



F I [g] = ^g \nl3(A^^ 



where 



In/?(A£,0) = lna+ [ — 

Jo y 



g(y) 
go 



(C13) 
H dyg(y) 

y go 

(C14) 



To evaluate integrals when A 1 <C |x| <C £, we intro- 
duce a new function 



h(u) 



<,g(u) J 1 as u — > oo, 

3oo l(go/.goo)w as u -> 0, 



(C18) 



where g^, is defined in Eq. l|3.57[l Then = 

ffoo (^/7r)(|x|/O X(x), where 



/■ A|x| rfw 

1= / ^-y-^/W-Uy)]. 

Jo y 



(C19) 



This integral has a potential infrared divergence as 
£/|x| — > oo when > 2. To isolate it, we break up 
the integral from to A|x| into one from to 1 and 
another from 1 to A|x|. There are no troubles with ul- 
traviolet divergences in the second integral, and in it, 
we can let A|x| — ► oo and replace h(y£/\x\) by its in- 
finite argument limit of one. In the integral from 
to 1, we extract the small y behavior of 1 — Jo(y) via 
1 - My) = (y 2 /4) + [1 - My) - (y 2 / 4 )]- The second 
part of this expression vanishes as y A at small y, and 
there is no infrared divergence in the integral involving 
it so long as <j> < 4. Thus, we have 

t = i x + f - My) - (y 2 /4)] (C20) 
Jo y 



dy_ 

i y 



y -+[\ - My)}, 



(C21) 



20 



where 



1 , £\<f>-2 r «/|x| 



4 VW 
4 x 



0-2 



C/|x| 



dyu 



duu 1 ^h(u) 
duu 1 -^>h(u) 



+ y u 1 -^^)-!]}. 



(C25) 



Using J-L 
(|3~(3l"|) with 



^2-0 _ j]/^ - </»), we arrive at Eq. 



dyy-V+Vil-Mv)] 
1 du. 



(C26) 



Lai/ = / —h{u)+ r — [h(u)-l] (C27) 
Jo u Ji u 



C 2 (4>) = / duu^hiu). 



(C28) 



The evaluation of the £ ^> |x| limit of F^fg] is straight- 
forward. The result is 



* Jo V 



(C29) 



where P n (x) is the nth-order Legendre Polynomial, we 





find 




(C22) 


T A = 










(C23) 


T B = 








7T 


(C24) 


T c = 





^[as] + -P2(cos6» x )/2[5s] 



— P 2 (cos 0* ) + - sin 2 ^ cos ' 



(C35) 
(C36) 

(C37) 
h[gc], (C38) 



where 



lob] = / <%(?£)[! -io(g|x|)] 



^[ff] = / d?fl(?0ia(?|x|) 



^[ff] = / (%(<7£)j 4 (g|x|). 
Jo 



Thus, 



+ J P 2 (C0S^)J 2 [52] 

— P 2 {cos8 x ) + - sin 2 X cos 4^ 



(C39) 
(C40) 
(C41) 



(C42) 
(C43) 

h\gc], (C44) 



where gi = g A + \gB - j^gc and g 2 = \gs + %y9C- Thus 
we need only evaluate the three integrals I±, I 2l and I3. 



2. Three Dimensions 

To evaluate the integrals T a in 3-D, we make use of the 
ZD plane-wave decomposition: 

00 1 
e^ x = 47r^^(g|x|) Y, Y lm (n q )YC m (n x ), (C30) 

1=0 m=-l 

where Q q = (9 q ,(j) q ) and £l x — (9 X ,4> X ) are, respectively, 
the polar angles of q and x, Yj m (f2) are spherical harmon- 
ics, and j n (u) is the nth order spherical Bessel function. 
Then, noting that 



q\ = sin 2 q 


= ^[l-P 2 (cos 


o q )], 


(C31) 


-2-2 -2 

Ixly = g sm 


(9,(1 -cos40 ? ) 




(C32) 




- 1OP 2 (cos0,) + 


3P 4 (cos 0,)] 


(C33) 


1 / 

~ 16 V 


^2 4 4![y 44 (fi ? )H 




(C34) 



a. |x| S> £ > A 1 in Three Dimensions 

In this limit, in integrals with integrands proportional 
to j n (q\x\), we set y = q\x\, replace g(y£/|x|) by g Q and 
replace the upper limit, A|x|, of integration by 00. In the 
part of the integral I± not proportional to j'o(g|x|), we set 
y = q£. The result is 



T 77 

h ~ 9o 2 



T 9o 
11 ~ 1 — r 



2 f H 9(M) dy _± 
k Jo 9o |x| 



(C45) 

(C46) 
(C47) 



r 3o f°° . , x 3.g 7r 
l x l Jo !6|x| 

b. A -1 <C |x| <C £ in Three Dimensions 



To treat this limit, as in 2D, we use the function h(u) 
[Eq. (EH]. To evaluate h, we break up the limits of 
integration in much the same way we did in 2D. The 
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result is 
h = .g co r N " 1 



(C48) 



J dyy-f[l-j Q {y)]- J dyy-*j (y) (C49) 



+ T — [(Alxl) 1 ^-!] 



(C50) 



1-0 r u 



duu 



\h{u)-l]) (C51) 



for < (j> < 3. The dominant behavior for 1 < <fi < 3 and 
< 6 < 1 is then 



where 



1< < 3 
< 1. 
(C52) 



M4>) = / 2T [1-Jo(2/)] (C53a) 



C 3 (0) = / y-^joiy). (C53b) 
Jo 



3. One Dimension 

In ID, there is only one function to evaluate 
Hx) = 21* j ^9(qm ~ cos(q\x\)] 

= — / dy-g(yt/\x\)(l-cosy). (C54) 
K Jo V 

The limit \x\ ^> £ is obtained as before by replacing 
g(y£/\x\) with g and letting A\x\ — > oo: 



'dy^— ^= ff o^W. (C55) 



In the limit A 1 <C |x| <C £, we introduce as in 
2D and 3D: T{x) = (2g oc \x\/n)]C, where 



K = 



dyh(y^/\x\)y-^(l-cosy) 



= /C x + / ^jT (2+ « [1 - cos y - (y 2 /2)] 



"(2+0) 



[1-cosy], (C56) 



where 



dyy-*h(yZ/\x\) 



2\\l\) {l-d>([\l\ 



(C57) 



duh(u)u 



-0 



duvT* 



[h(u)-l]}. 



Combining Eqs. (|C56|I with (|C57fl . we find 

f^ockr+Mi^) if < i, 

F(x)~l±liut/\x\ if 0=1, (C58) 

U^'M 2 ^) if > l, 

where z/ is given by Eq. ( IC27() and C2 ((f)) is given by Eq. 
(IC28II and where 



1 — cos y 



y 



,2+0 



(C59) 



1 ~ \ 1 ; A? d W5 H - soor'M-^csfo) < < 1. ' 



The |x| <C £ limits of both 22 and -Z3 can be obtained 
by simply by replacing g(q£) by {qO^doo- 

h ~ gooC*]^- 1 dyy-^j 2 (y) < < 3 (C60) 



h-gooC^- 1 dyy-*u(y) < < 5. (C61) 
Jo 



APPENDIX D: EVALUATION OF &,(x) 

In this appendix we will evaluate the rotational corre- 
lation function G w (x) in two dimensions. To lowest order 
in A K , 

G w (q) = \tri£rii>qrq r ix%X<i)xl 1 A-Q L )^2s app >, (Dl) 



where S app ' is defined in Eq. (|B4J| . The product e r i£ r 'i' is 
simply S rr 'Sur — 6 r ir5i r >, and e r ie r 'i'q r q' r — q 2 5j i ,. When 
this operates on Xi P Xi> P >> it projects out the transverse 
part leaving Sp p ,/(fi 2 q 2 ). Thus 

G -(q) = %-J2 S * = -%( A i4 - (D2) 

where = A 3 + 4A 6 + 2A 7 + 2A 8 and A£ = 4A 3 + 
15A 6 + 4A 7 . 

When A^ c (q) have a Lorentizan form, we need to 
evaluate two integrals to determine G w (x) : 



and 

*a(x) = 



A 9x^_ iq . x 



(2tt) 2 g 2 + 



2W0 g 2 + r 2 
1 ko(|x|/0 



2tt 



■ cos 2 sin 2 e *9|x| cos(0-v.) 



16tt 



cos 4^ 



48£ 4 , 4£ 



+ T ^+jr 4 (|x|/0 , (D4) 



where q = g(cos (j>, sin </>) , x = |x| (cos tp, sin ■)/>) and K n (y) 
is the Bessel function of imaginary argument. 
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